%%bash
wget http://kodomo.cmm.msu.ru/~golovin/bilayer/dppc.itp
wget http://kodomo.cmm.msu.ru/~golovin/bilayer/lipid.itp
wget http://kodomo.cmm.msu.ru/~golovin/bilayer/dppc.gro
wget http://kodomo.cmm.msu.ru/~golovin/bilayer/b.top
wget http://kodomo.cmm.msu.ru/~golovin/bilayer/em.mdp
wget http://kodomo.cmm.msu.ru/~golovin/bilayer/pr.mdp
wget http://kodomo.cmm.msu.ru/~golovin/bilayer/md.mdp
! genconf -f dppc.gro -o b_64.gro -nbox 4 4 4
! editconf -f dppc.gro -o dppc.pdb
! editconf -f b_64.gro -o b_64.pdb
Посмотрим в PyMol
Посмотри 64 молекулы
! editconf -f b_64.gro -o b_ec -d 0.5
! grompp -f em -c b_ec -p b -o b_em -maxwarn 2
! mdrun -deffnm b_em -v
Посмотрим изменение силы
Fmax=4.37970e+05
Fmax=6.16887e+02
! genbox -cp b_em -p b -cs spc216 -o b_s
! grompp -f pr -c b_s -p b -o b_pr -maxwarn 1
! mdrun -deffnm b_pr -v
! grompp -f em -c b_s -p b -o b_empr -maxwarn 1
! mdrun -deffnm b_empr -v
! grompp -f pr -c b_empr -p b -o b_pr -maxwarn 1
!mdrun -deffnm b_pr -v
Переформатируем в pdb
! editconf -f b_pr.gro -o b_pr.pdb
! editconf -f b_s.gro -o b_s.pdb
Посмотрим b_pr
И теперь b_s
В целом, существенной разницы нет, возможно, есть небольшая разница в количество и положениях молекул воды.
Теперь посмотрим результат моделирования
! echo 2 | trjconv -f b_md.xtc -s b_md.tpr -o b_pbc_1.pdb -skip 20
Попробуем второй способ визуализации
! echo 2 | trjconv -f b_md.xtc -s b_md.tpr -o b_pbc_1.pdb -skip 20 -pbc mol
Посмотрим, что получилось
from IPython.display import Video
Video('b_pbc_1.mpg')
Видим, что довольно быстро, уже в районе 20 кадра образуется бислой.
! echo 0 | g_traj -f b_md.xtc -s b_md.tpr -ob box_1.xvg
Прочитаем результат и построим результат с учётом того, что нормалью к поверхности является ось X
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
%matplotlib inline
data = np.loadtxt('box_1.xvg', comments=['#', '@'])[:,:4]
plt.plot(data[:, 0], data[:, 1] * data[:, 3] / 64)
plt.title('square')
plt.xlabel('time')
Text(0.5,0,u'time')
Как видно по графику, нормализованная площадь падает со временем, это значит, что действительно происходит сборка липидного слоя.
Теперь построим график изменения гидрофобной и гидрофильной поверхности
! g_sas -f b_md.xtc -s b_md.tpr -o sas_b.xvg -dt 100
data = np.loadtxt('sas_b.xvg', comments=['#', '@'])
time = data[:,0]
plt.plot(time, data[:,1], color='green')
plt.plot(time, data[:,2], color='blue')
patch1 = mpatches.Patch(color='green', label='Hydrophobic')
patch2 = mpatches.Patch(color='blue', label='Hydrophilic')
plt.legend(handles=[patch1,patch2])
plt.title('Hydrophobic and hydrophilic surface square')
plt.xlabel('Time')
Text(0.5,0,u'Time')
Видим, что площадь гидрофобных поверхностей падает быстрее, чем площадь гидрофильных, это ожидаемо, так как в процессе сборки липидного слоя молекулы поворачиваются гидрофобными частями внутрь слоя, так как контакты гидробочных поверхностей с молекулами воды энергетически невыгодны. Площадь гидрофильных поверхностей тоже немного падает за счёт компактификации струкруры.
Теперь посмотрим на изменение меры порядка.
!wget http://kodomo.cmm.msu.ru/~golovin/bilayer/sn1.ndx
! g_order -s b_md -f b_md.xtc -o ord_end.xvg -n sn1.ndx -b 45000 -d X
! g_order -s b_md -f b_md.xtc -o ord_start.xvg -n sn1.ndx -e 5000 -d X
import pandas as pd
start = pd.read_csv('ord_start.xvg', header = None, skiprows = 12, sep='\s+')
end = pd.read_csv('ord_end.xvg', header = None, skiprows = 12, sep='\s+')
plt.plot(start[0], start[3], marker='.', label='start')
plt.plot(end[0], end[3], marker='.', label='end')
plt.xlabel('atom')
plt.ylabel('order')
plt.legend()
<matplotlib.legend.Legend at 0x7fd1906ad910>
Видим, что порядок в результате сборки в конце увеличилась мера порядка, это значит, что в результате сборки структура стала более упорядоченной.